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Abstract 



Generalized Linear Models (GLMs) and Single Index Models (SIMs) provide powerful generalizations 
of linear regression, where the target variable is assumed to be a (possibly unknown) 1-dimensional 
function of a linear predictor. In general, these problems entail non-convex estimation procedures, and, 
in practice, iterative local search heuristics are often used. Kalai and Sastry (2009) recently provided 
the first provably efficient method for learning SIMs and GLMs, under the assumptions that the data 
are in fact generated under a GLM and under certain monotonicity and Lipschitz constraints. However, 
to obtain provable performance, the method requires a fresh sample every iteration. In this paper, 
we provide algorithms for learning GLMs and SIMs, which are both computationally and statistically 
efficient. We also provide an empirical study, demonstrating their feasibility in practice. 

1 Introduction 

The oft used linear regression paradigm models a target variable Y as a linear function of a vector-valued 
input X. Namely, for some vector w, we assume that E[F|X] = w ■ X. Generalized linear models (GLMs) 
provide a flexible extension of linear regression, by assuming the existence of a "link" function g such 
that E[y|X] = g^ 1 (w ■ X). g "links" the conditional expectation of Y to X in a linear manner, i.e. 
g(K[Y\X]) — w ■ X (see [MN89 for a review). This simple assumption immediately leads to many practical 
models, including logistic regression, the workhorse for binary probabilistic modeling. 

Typically, the link function is assumed to be known (often chosen based on problem-specific constraints) , 
and the parameter w is estimated using some iterative procedure. Even in the setting where g is known, 
we are not aware of a classical estimation procedure which is computationally efficient, yet achieves a good 
statistical rate with provable guarantees. The standard procedure is iteratively reweighted least squares, 
based on Newton- Ralphson (see |MN89j ). 

In Single Index Models (SIMs), both g and w are unknown. Here, we face the more challenging (and 
practically relevant) question of jointly estimating g and w, where g may come from a large non-parametric 
family such as all monotonic functions. There are two issues here: 1) What statistical rate is achievable 
for simultaneous estimation of g and w? 2) Is there a computationally efficient algorithm for this joint 
estimation? With regards to the former, under mild Lipschitz-continuity restrictions on <? _1 , it is possible to 
characterize the effectiveness of an (appropriately constrained) joint empirical risk minimization procedure. 

*This work was partially supported by grants NSF-CCF-04-27129 and NSF-CCF-09-64401. 



1 



This suggests that, from a purely statistical viewpoint, it may be worthwhile to attempt to jointly optimize 
g and w on the empirical data. 

However, the issue of computationally efficiently estimating both g and w (and still achieving a good 
statistical rate) is more delicate, and is the focus of this work. We note that this is not a trivial problem: 
in general, the joint estimation problem is highly non-convex, and despite a significant body of literature 
on the problem, existing methods arc usually based on heuristics, which are not guaranteed to converge 
to a glob al optimum (see for instance |WHI93l IHH941 IMHS981 INT041 IRWY08] ). We note that recently, 
SSSSlD] presented a kernel-based method which does allow (improper) learning of certain types of GLM's 
and SIM's, even in an agnostic setting where no assumptions are made on the underlying distribution. On 
the flip side, the formal computational complexity guarantee degrades super-polynomially with the norm of 
w, which SSSS10J show is provably unavoidable in their setting. 

The recently proposed Isotron algorithm [KS09 provides the first provably efficient method for learning 
GLMs and SIMs, under the common assumption that g^ 1 is monotonic and Lipschitz, and assuming the data 
corresponds to the model. The algorithm attained both polynomial sample and computational complexity, 
with a sample size dependence that does not depend explicitly on the dimension. The algorithm is a variant 
of the "gradient-like" perceptron algorithm, with the added twist that on each update, an isotonic regression 
procedure is performed on the linear predictions. Recall that isotonic regression is a procedure which finds 
the best monotonic one dimensional regression function. Here, the well-known Pool Adjacent Violator (PAV) 
algorithm provides a computationally efficient method for this task. 

Unfortunately, a cursory inspection of the Isotron algorithm suggests that, while it is computationally 
efficient, it is very wasteful statistically, as each iteration of the algorithm throws away all previous training 
data and requests new examples. Our intuition is that the underlying technical reasons for this are due to 
the fact that the PAV algorithm need not return a function with a bounded Lipschitz constant. Furthermore, 
empirically, it not clear how deleterious this issue may be. 

This work seeks to address these issues both theoretically and practically. We present two algorithms, the 
GLM-tron algorithm for learning GLMs with a known monotonic and Lipschitz <? _1 , and the L-lsotron algo- 
rithm for the more general problem of learning SIMs, with an unknown monotonic and Lipschitz g . Both 
algorithms are practical, parameter-free and are provably efficient, both statistically and computationally. 
Moreover, they are both easily kernelizable. In addition, we investigate both algorithms empirically, and 
show they are both feasible approaches. Furthermore, our results show that the original Isotron algorithm 
(ran on the same data each time) is perhaps also effective in several cases, even though the PAV algorithm 
does not have a Lipschitz constraint. 

More generally, it is interesting to note how the statistical assumption that the data are in fact generated 
by some GLM leads to an efficient estimation procedure, despite it being a non-convex problem. Without 
making any assumptions, i.e. in the agnostic setting, this problem is at least hard as learning parities with 
noise. 

2 Setting 

We assume the data (x,y) are sampled i.i.d. from a distribution supported on B^ x [0, 1], where B^ = {x £ 
R d : \\x\\ < 1} is the unit ball in c?-dimensional Euclidean space. Our algorithms and analysis also apply 
to the case where B^ is the unit ball in some high (or infinite)-dimensional kernel feature space. We assume 
there is a fixed vector w, such that \\w\\ < W, and a non-decreasing 1-Lipschitz function u : R — > [0,1], such 
that E[y|&] = u{w-x) for all x. Note that u plays the same role here as g^ 1 in generalized linear models, and 
we use this notation for convenience. Also, the restriction that u is 1-Lipschitz is without loss of generality, 
since the norm of w is arbitrary (an equivalent restriction is that ||u>|| = 1 and that u is VF-Lipschitz for an 
arbitrary W). 

Our focus is on approximating the regression function well, as measured by the squared loss. For a real 
valued function h : — > [0, 1], define 

err(/i) = E {x>y) [(h(x) - y) 2 } 
e(h) = err(h) — err(E[y\x\) 

= E( x ,y) [(h(x) - u(w ■ x)) 2 ] 
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Algorithm 1 GLM-tron 



Input: data ((i i5 1 e R d x [0, 1], tt : M -> [0, 1]. 
:= 0; 

for * = 1,2,... do 
h*(x) := «(«>* • a;); 

u> t+1 := to* H Y^G/i - "(to* • Xi))xi\ 

m *■ — ' 

i—l 

end for 



err(/i) measures the error of h, and s(h) measures the excess error of h compared to the Bayes-optimal 
predictor x H- u(w ■ x). Our goal is to find h such that e(h) (equivalently, err(/i)) is as small as possible. 
In addition, we define the empirical counterpart efr(/i), i(h), based on a sample (x%, yi), . . . , (x m , y m ), to 

be 

^ m 

erf(/i) = — V(M^) - y,) 2 
rn * — ' 

8=1 

e(/i) = — VY/ife) - u(tu • Xi)) 2 . 

i=l 

Note that e is the standard fixed design error (as this error conditions on the observed x's). 

Our algorithms work by iteratively constructing hypotheses h* of the form /i'(x) = ^(w* ■ x), where u* 
is a non-decreasing, 1-Lipschitz function, and w* is a linear predictor. The algorithmic analysis provides 
conditions under which is small, and using statistical arguments, one can guarantee that e{h l ) would 

be small as well. 

To simplify the presentation of our results, we use the standard O(-) notation, which always hides only 
universal constants. 



3 The GLM-tron Algorithm 

We begin with the simpler case, where the transfer function u is assumed to be known (e.g. a sigmoid), 
and the problem is estimating w properly. We present a simple, parameter-free, perceptron-like algorithm, 
GLM-tron, which efficiently finds a close-to-optimal predictor. We note that the algorithm works for arbitrary 
non-decreasing, Lipschitz functions u, and thus covers most generalized linear models. The pseudo-code 
appears as Algorithm [T] 

To analyze the performance of the algorithm, we show that if we run the algorithm for sufficiently 
many iterations, one of the predictors h obtained must be nearly-optimal, compared to the Bayes-optimal 
predictor. 

Theorem 1. Suppose (x\, yi), . . . , (x m , y m ) are drawn independently from a distribution supported on x 
[0, 1], such that E[j/|x] = u(w ■ x), where \\w\\ < W , and u : K. — > [0, 1] is a known non- decreasing 1-Lipschitz 
function. Then for any 8 € (0, 1), the following holds with probability at least 1 — 5: there exists some 
iteration t < 0(W \Jmj log(l/5)) of GLM-tron such that the hypothesis h l {x) = u(w l ■ x) satisfies 

mBx{t{h%e{h*)} < O . 

In particular, the theorem implies that some h l has e(/i*) = 0(l/^/m). Since e(h f ) equals err(/i,') up to 
a constant, we can easily find an appropriate h l by using a hold-out set to estimate err(/r), and picking the 
one with the lowest value. 

The proof is along similar lines (but somewhat simpler) than the proof of our subsequent Thm. [2] 
The rough idea of the proof is showing that at each iteration, if s(h ) is not small, then the squared 
distance — u>*|| is substantially smaller than — w\\ 2 . Since this is bounded below by 0, and 
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Algorithm 2 L-lsotron 



Input: data ((x^y^))^ G R d x [0, 1]. 
w 1 :=0; 

for * = 1,2,... do 

u f := LPAV((w' • xi,yi), . . . , (w 1 ■ x m ,y m )) 

^ m 

w t+1 := w t + — y^(yi - u t (w t ■ x i ))x l 

i=l 

end for 



w° — to 1 1 < W 2 , there is an iteration (arrived at within reasonable time) such that the hypothesis h* at 
that iteration is highly accurate. The proof is provided in Appendix | A. 1| 

4 The L-lsotron Algorithm 

We now present L-lsotron, in Algorithm|2j which is applicable to the harder setting where the transfer function 
u is unknown, except for it being non-decreasing and 1-Lipschitz. This corresponds to the semi-parametric 
setting of single index models. 

The algorithm that we present is again simple and parameter-free. The main difference compared to 
GLM-tron algorithm is that now the transfer function must also be learned, and the algorithm keeps track 
of a transfer function u which changes from iteration to iteration. The algorithm is also rather similar to 
the Isotron algorithm |KS09j . with the main difference being that instead of applying the PAV procedure 
to fit an arbitrary monotonic function at each iteration, we use a different procedure, LPAV, which fits a 
Lipschitz monotonic function. This difference is the key which allows us to make the algorithm practical 
while maintaining non-trivial guarantees (getting similar guarantees for the Isotron required a fresh training 
sample at each iteration). 

The LPAV procedure takes as input a set of points (z\, yi), . . . , (z m , y m ) in M 2 , and fits a non-decreasing, 
1-Lipschitz function u, which minimizes X)i=i( u ( z i) — Vi) 2 ■ This problem has been studied in the literature, 
and we followed the method of [YW09] in our empirical studies. The running time of the method proposed 
in |YW09| is 0(m 2 ). While this can be slow for large-scale datasets, we remind the reader that this is a 
one-dimensional fitting problem, and thus a highly accurate fit can be achieved by randomly subsampling 
the data (the details of this argument, while straightforward, are beyond the scope of the paper). 

We now turn to the formal analysis of the algorithm. The formal guarantees parallel those of the 
previous subsection. However, the rates achieved are somewhat worse, due to the additional difficulty of 
simultaneously estimating both u and w. It is plausible that these rates are sharp for information-theoretic 
reasons, based on the 1-dimensional lower bounds in jZha02aj (although the assumptions are slightly different, 
and thus they do not directly apply to our setting). 

Theorem 2. Suppose {x\, yi), . . . , (x m , y m ) are drawn independently from a distribution supported on x 
[0,1], such that E[y|x] = u{w ■ x), where \\w\\ < W, and u : K — > [0,1] is an unknown non- decreasing 
1-Lipschitz function. Then the following two bounds hold: 

1. (Dimension- dependent) With probability at least 1—5, there exists some iteration t < O 
of L-lsotron such that 

max{e(/i t ), eQv*)} <°[[ — — J 

2. (Dimension-independent) With probability at least 1—8, there exists some iteration t < O 
of L-lsotron such that 

max { ^), £ (^) } <0^ 2l0 ^) 1/4 ) 



/ Wm \ 
\d\og(Wm/8) ) 



1/3 



U log(m/5)J J 



4 



As in the case of Thm. [T] one can easily find h l which satisfies the theorem's conditions, by running the 
L-lsotron algorithm for sufficiently many iterations, and choosing the hypothesis h which minimizes err(/i') 
based on a hold-out set. 



5 Proofs 

5.1 Proof of Thm. H 

First we need a property of the LPAV algorithm that is used to find the best one-dimensional non-decreasing 
1-Lipschitz function. Formally, this problem can be defined as follows: Given as input ({z,, Ui})"^ € 
[— W, W] x [0, 1] the goal is to find y%, . . . ,y m such that 



-X>-*) 2 ' (!) 



m 
»=i 

is minimal, under the constraint that yi — u(zi) for some non-decreasing 1-Lipschitz function u : [— W, W] >-> 
[0,1]. After finding such values, LPAV obtains an entire function u by interpolating linearly between the 
points. Assuming that Zi are in sorted order, this can be formulated as a quadratic problem with the following 
constraints: 

m - Vi+x < 1 < i < m (2) 

Vi+i ~y%~ (Zi+i - Zi) <0 1 < i < m (3) 

Lemma 1. Let (zx, yi), . . . , (z m , y m ) be input to LPAV where Zi are increasing and yi € [0, 1]. Let y\, . . . , y m 
be the output of LPAV. Let f be any function such that /(/?) — f{a) > ft — a, for j3 > a, then 

m 
i=l 

Proof. We first note that (yj —yj) — 0, since otherwise we could have found other values for yi, . . . ,y m 

which make ([T]) even smaller. So for notational convenience, let yo = 0, and we may assume w.l.o.g. that 
f(Vo) = 0- Define a t = Y^iiVj ~ Vj)- Then we have 



^2ivi - y%)(f(vi) - Zi) = 

rn 

E^(Cf(w)-««)-(/(w-x)-«<-x))- (4) 



i=i 



Suppose that oi < 0. Intuitively, this means that if we could have decreased all values . . . y m by an 
infinitesimal constant, then the objective function ([I]) would have been reduced, contradicting the optimality 
of the values. This means that the constraint y~i — yi+\ < must be tight, so we have (f(yi+i) — z i+ i) — 
(fiVi) ~ z i) = — z i+i +Zi < (this argument is informal, but can be easily formalized using KKT conditions). 
Similarly, when o~i > 0, then the constraint jji+i — yi~- {zi+i — Zi) < must be tight, hence /(?/;+i) — /($«) > 
Vi+i ~ Hi — ( z i+i ~ z i) —\ 0. So in either case, each summand in Q must be non-negative, leading to the 
required result. □ 

We also use another result, for which we require a bit of additional notation. At each iteration of the 
L-lsotron algorithm, we run the LPAV procedure based on the training sample (xi,yi), . . . , (x m , y m ) and the 
current direction w , and get a non-decreasing Lipschitz function u . Define 



Vi y\ = 



x j 



Recall that w, u are such that E[y|x] = u(w-x), and the input to the L-lsotron algorithm is (xi, y±), . . . , (x TO , y m )- 
Define 

Vi yi = u(w ■ Xi) 
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to be the expected value of each yi. Clearly, we do not have access to yi. However, consider a hypothetical 
call to LPAV with inputs ((w* ■ £i, yi))fc=i> and suppose LPAV returns the function u*. In that case, define 



Vi y\ = • x{). 

for all i. Our proof uses the following proposition, which relates the values y\ (the values we can actually 
compute) and y\ (the values we could compute if we had the conditional means of each j/j). The proof of 
Proposition [T] is somewhat lengthy and requires additional technical machinery, and is therefore relegated to 
Appendix [B| 

Proposition 1. With probability at least 1 — 5 over the sample {{xi,yi)}™ =1 , it holds for any t that 
m T^iLi \ vl ~ Vi\ * s a t m °st the minimum of 



O 



dW 2 \og(Wm/5) 



1/3N 



and 



O 



W 2 \og{m/5) 



1/4N 



The third auxiliary result we'll need is the following, which is well-known (see for example |STC04| . 
Section 4.1). 



Lemma 2. Suppose Z\,...,z m are i.i.d. 0-mean random variables in a Hilbert space, such that Pr(||:rj|| < 
1) = 1. Then with probability at least 1 — 5, 



m ^ 



< 2 



1 



v/log(lA5)/2 \ 



With these auxiliary results in hand, we can now turn to prove Thm. [2] itself. The heart of the proof is 
the following lemma, which shows that the squared distance \\w — w\\ between w and the true direction 
w decreases at each iteration at a rate which depends on the error of the hypothesis e(/i*): 



Lemma 3. Suppose that -w\\ <W and \\{l/m) Y^LxiVi ~ Vi) x i\\ < Vi and (1/m) YaLx \vl ~ vl\ <m- 
Then 



I t II 2 
\w — w\\ 



■w\\ >e{h t )-bW{r il +ri2) 



Proof. We have 



IV 



i+i 



w + w — w\ 



= \\w t+1 - II + ||w* - w\\ + 2{w t+1 - w 1 ) ■ (w 1 - w) 



Since w t+ — u>* = {1/m) 53i=i(fi — y\) x ii substituting this above and rearranging the terms we get, 



I t II 2 
\w — w\\ 



\w t+1 — w\ 



11} 



^2(Vi - Vi)( w ■ %i ~ wt ■ x i) - 



^ m 

— yZiVi - y\)xi 



(5) 



G 



Consider the first term above, 

2 



to . 
i— 1 



^2(yi - yt)( w ■x l -w t ■ Xi) 



2 



TO 



X^ - ^) 21 * yiw-w*) (6) 



2 



TO 



2 



-X^-tfXwsi-w'-a*) (8) 



TO 



The term |6j is at least —2Wr\\, the term (|8| is at least —2Wi]2 (since — w*) < W). We thus consider 
the remaining term Q . Letting u be the true transfer function, suppose for a minute it is strictly increasing, 
so its inverse u" 1 is well defined. Then we have 

2 m 

— 22(yi ~ yt)( w -Xi-w* ■ Xi) 

i=l 

2 rn 



TO 

i=l 



„ m 
m ^ — ' 



TO 

i=l 



The second term in the expression above is positive by Lemma [T] As to the first term, it is equal to 
I; Ei=i(yi-yf)(' u ~ 1 fe)-' lt " 1 (y|)) J which by the Lipschitz property of u is at least £ E™ ife-2/') 2 = 2e(ji*). 
Plugging this in the above, we get 



- - #)(«> ■ Xi - w* ■ x t ) > 2e(h t ) - 2W{ m + m ) 

111 ^— » 



(9) 



j=i 



This inequality was obtained under the assumption that u is strictly increasing, but it is not hard to verify 
that the same would hold even if u is only non-decreasing. 

The second term in ([5| can be bounded, using some tedious technical manipulations (see (14) and (15) 
in the supplementary material), by 



TO ^-^ 

Combining ^ and (10)) in <f5j> , we get 



i=l 



I t II 2 



„t+i 



w\\ > 2i(h t ) - - VK(5r?i + 2r? 2 ) 



(10) 



(11) 



Now, we claim that 



ift) - i(h*) > £|$-$|>-2%, 
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since 



1 m 

i=l 
1 m 

= ^£(fff-tf + fij-fc) a 

i=l 
1 m 

= i£(fif-ft) 9 

i=l 

V i=l / 



and we have that \y\ + y\ — 2y,| < 2. Plugging this into ( |ll[ ) leads to the desired result. □ 

The bound on in Thm. [2]now follows from Lemma[3] Using the notation from Lemma[3j r\\ can be 
set to the bound in Lemma|2j since {(y$ — yi)xi}^ 1 are i.i.d. 0-mean random variables with norm bounded 

by 1. Also, ri2 can be set to any of the bounds in Proposition fT] 772 is clearly the dominant term. Thus, we 

Si 1 2 2 — 

holds, so either \\w t+1 - w\\ < Ww* - w\\ - W(r)i + 772), or e(/V) < 3VF(r7i + 772). If the 

latter is the case, we are done. If not, since ||u> t+1 — toll > 0, and \\w° — w\\ 2 — \\w\\ 2 < W 2 , there can be 
at most W 2 /(W(i]i + i] 2 )) = W/(r]i + rj 2 ) iterations before i(h l ) < 6Wi]. Plugging in the values for 771,772 
results in the bound on 

Finally, to get a bound on s(h l ), we utilize the following uniform convergence lemma: 

Lemma 4. Suppose thatE[y\x] — u((w,x)) for some non- decreasing 1-Lipschitzu andw such that \\w\\ < W. 
Then with probability at least 1 — 5 over a sample (x±,yi), . . . , (x m ,y m ), the following holds simultaneously 
for any function h{x) — u(w ■ x) such that \\w\\ < W and a non- decreasing and 1-Lipschitz function u: 



\e(h)-i(h)\ < O 



'W 2 \og{m/5) 



The proof of the lemma uses a covering number argument, and is shown as part of the more general 
Lemma [7] in the supplementary material. This lemma applies in particular to h l . Combining this with the 
bound on £(/i 4 ), and using a union bound, we get the result on s(h ) as well. 



6 Experiments 

In this section, we present an empirical study of the GLM-tron and the L-lsotron algorithms. The first 
experiment we performed is a synthetic one, and is meant to highlight the difference between L-lsotron and 
the Isotron algorithm of [KS09J. In particular, we show that attempting to fit the transfer function without 
any Lipschitz constraints may cause Isotron to overfit, complementing our theoretical findings. The second 
set of experiments is a comparison between GLM-tron, L-lsotron and several competing approaches. The 
goal of these experiments is to show that our algorithms perform well on real-world data, even when the 
distributional assumption required for their theoretical guarantees does not precisely hold. 

6.1 L-lsotron vs Isotron 

As discussed earlier, our L-lsotron algorithm (Algorithm [2]) is similar to the Isotron algorithm of }KS09| . with 
two main differences: First, we apply LPAV at each iteration to find the best Lipschitz monotonic function 
to the data, while they apply the PAV (Pool Adjacent Violator) procedure to fit a monotonic (generally 
non-Lipschitz) function. The second difference is the theoretical guarantees, which in the case of Isotron 
required working with a fresh training sample at each iteration. 

While the first difference is inherent, the second difference is just an outcome of the analysis. In particular, 
one might still try and apply the Isotron algorithm, using the same training sample at each iteration. While 
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Figure 1: The link function as predicted by LIsotron (blue) and Isotron (red). The domain of both functions 
was normalized to [—1,1]. 



we do not have theoretical guarantees for this algorithm, it is computationally efficient, and one might 
wonder how well it performs in practice. As we see later on, it actually performs quite well on the datasets 
we examined. However, in this subsection we provide a simple example, which shows that sometimes, the 
repeated fitting of a nora-Lipschitz function, as done in the Isotron algorithm, can cause overfit and thus 
hurt performance, compared to fitting a Lipschitz function as done in the L-lsotron algorithm. 

We constructed a synthetic dataset as follows: In a high dimensional space (d = 400), we let w = 
(1, 0, . . . , 0) be the true direction. The transfer function is u(t) = (l+<)/2. Each data point x is constructed as 
follows: the first coordinate is chosen uniformly from the set {—1, 0, 1}, and out of the remaining coordinates, 
one is chosen uniformly at random and is set to 1. All other coordinates are set to 0. The y values are chosen 
at random from {0, 1}, so that E[y|a;] = u(w ■ x). We used a sample of size 600 to evaluate the performance 
of the algorithms. 

In the synthetic example we construct, the first attribute is the only relevant attribute. However, because 
of the random noise in the y values, Isotron tends to overfit using the irrelevant attributes. At data points 
where the true mean value u(w-x) equals 0.5, Isotron (which uses PAV) tries to fit the value or 1, whichever 
is observed. On the other hand, L-lsotron (which uses LPAV) predicts this correctly as close to 0.5, because of 
the Lipschitz constraint. Figure [T] shows the link functions predicted by L-lsotron and Isotron on this dataset. 
Repeating the experiment 10 times, the error of L-lsotron, normalized by the variance of the y values, was 
0.338±0.058, while the normalized error for the Isotron algorithm was 0.526±0.175. In addition, we observed 
that L-lsotron performed better rather consistently across the folds - the difference between the normalized 
error of Isotron and L-lsotron was 0.189 ± 0.139. 



Table 1: Mean squared error normalized by the variance (mean and standard deviation across 10 folds). 



dataset L-Iso GLM-t Iso Lin-R Log-R SIM 



communities 


0.34 


± 


0.04 


0.34 


± 


0.03 


0.35 


± 


0.04 


0.35 


± 


0.04 


0.34 


± 


0.03 


0.36 


± 


0.05 


concrete 


0.35 


± 


0.06 


0.40 


± 


0.07 


0.36 


± 


0.06 


0.39 


± 


0.08 


0.39 


± 


0.08 


0.35 


± 


0.06 


housing 


0.27 


± 


0.12 


0.28 


± 


0.11 


0.27 


± 


0.12 


0.28 


± 


0.12 


0.27 


± 


0.11 


0.26 


± 


0.09 


parkinsons 


0.89 


± 


0.04 


0.92 


± 


0.04 


0.89 


± 


0.04 


0.90 


± 


0.04 


0.90 


± 


0.04 


0.92 


± 


0.03 


winequality 


0.78 


± 


0.07 


0.81 


± 


0.07 


0.78 


± 


0.07 


0.73 


± 


0.08 


0.73 


± 


0.08 


0.79 


± 


0.07 



6.2 Real World Datasets 

We now turn to describe the results of experiments performed on several UCI datasets. We chose the 
following 5 datasets: communities, concrete, housing, parkinsons, and wine-quality. 

On each dataset, we compared the performance of L-lsotron (L-Iso) and GLM-tron (GLM-t) with Isotron 
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Table 2: Performance comparison of L-lsotron with the other algorithms. The values reported are the 
difference in the normalized squared errors (mean and standard deviation across the 10 folds). Negative 
values indicate better performance than L-lsotron. 



dataset GLM-t Iso Lin-R Log-R SIM 



communities 


0.01 


± 


0.01 


0.01 


± 


0.02 


0.02 


± 


0.02 


0.01 


± 


0.02 


0.02 


± 


0.03 


concrete 


0.04 


± 


0.03 


0.00 


± 


0.01 


0.04 


± 


0.03 


0.04 


± 


0.03 


0.00 


± 


0.02 


housing 


0.02 


± 


0.05 


0.00 


± 


0.07 


0.02 


± 


0.05 


0.01 


± 


0.05 


-0.01 


± 


0.06 


parkinsons 


0.03 


± 


0.02 


0.00 


± 


0.01 


0.02 


± 


0.01 


0.02 


± 


0.01 


0.04 


± 


0.04 


winequality 


0.03 


± 


0.02 


0.00 


± 


0.01 


-0.05 


± 


0.03 


-0.05 


± 


0.03 


0.01 


± 


0.01 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



(a) concrete (b) communities 

Figure 2: The transfer function u as predicted by L-lsotron (blue) and Isotron (red) for the concrete and 
communities datasets. The domain of both functions was normalized to [—1, 1]. 

and several other algorithms. These include standard logistic regression (Log-R), linear regression (Lin-R) 
and a simple heuristic algorithm (SIM) for single index models, along the lines of standard iterative maximum- 
likelihood procedures for these types of problems (e.g., |Cos83j ). The algorithm works by iteratively fixing 
the direction w and finding the best transfer function u, and then fixing u and optimizing w via gradient 
descent. For each of the algorithms we performed 10-fold cross validation, using 1 fold each time as the test 
set, and we report averaged results across the folds. 

Table [I] shows the mean squared error of all the algorithms across ten folds normalized by the variance 
in the y values. Table [2] shows the difference between squared errors between the algorithms across the 
folds. The results indicate that the performance of L-lsotron and GLM-tron (and even Isotron) is comparable 
to other regression techniques and in many cases also slightly better. This suggests that these algorithms 
should work well in practice, while enjoying non-trivial theoretical guarantees. 

It is also illustrative to see how the transfer functions found by the two algorithms, L-lsotron and Isotron, 
compare to each other. In Figure [2j we plot the transfer function for concrete and communities. The 
plots illustrate the fact that Isotron repeatedly fits a non-Lipschitz function resulting in a piecewise con- 
stant function, which is less intuitive than the smoother, Lipschitz transfer function found by the L-lsotron 
algorithm. 
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A Appendix 

A.l Proof of Thm. Q] 

The reader is referred to GLM-tron (Alg. [I]) for notation used in this section. 

The main lemma shows that as long as the error of the current hypothesis is large the distance of our 
predicted direction vector w from the ideal direction w decreases. 

Lemma 5. At iteration t in GLM-tron, suppose || w* — w\\ < W, then if \\(l/m) Y)^^ (Vj — u(w ■ Xi))xi\\ w < 
r\, then 



„t+i 



w\\ > e{h l ) - 5Wn 



Proof. We have 



\w — w\ 



w t+1 - w\ 



y^XVi ~ U ( wt ' X i))( W ■ X l -W t ■ Xi) 



1=1 



^ m 

~~ YliVi -11(111* -Xi))o 



i=l 



(12) 



Consider the first term above 
2 



^2(y l -u(w t -x i ))( 



W-Xi — W t -Xi 



) = — y^(u(w-Xi)-u(w t -Xi))(w-x i -w t -Xi)-\ f - u(w ■ xA)Xi ) ■(tu- 

rn L — ' m \ f — ' / 



i=i 



Using the fact that u is non-decreasing and 1-Lipschitz (for the first term) and \w — < W and 
|| (l/rn) Y]a—i (Vj — u(w ■ Xi))xi\\ < ij, we can lower bound this by 



— ^(u(w-x l )-u(w t -Xi)) 2 -2Wn > 2i(h t ) - 2W<n. 

711 — ^ 



(13) 



i=l 



For the second term in (12 1, we have 

2 



rn 



< 



E(Vi ~ u(w t ■ x{))xi 

= 1 

1 rn 

E(yi - u(w • Xi)) 

—1 

^ m 

E(u(w • Xi) - u{w l • Xi)) 



1 m 

E(yi - u(w ■ x^ + u(w ■ Xi) - • Xi))a 



1 m ^ m 

E(l/i - u(w ■ X i ))x l X — y^(u(w ■ Xi) - U(w f ■ Xi))Xi 
m f * 



(14) 



Using the fact that ||(l/m) — u ( w ' 2P»))a?* || _• Vi an d using Jensen's inequality to show that 

||(l/m) Y%Li( u ( w ' x i) ~ u ( wt ' x i)) x i\f < (1/™) Y^iLi( u ( w ' x i) ~ u ( wt ' x i)) 2 = an d assuming W > 1, 

we get 



^ m 

— ^(Vi - U ( W ■ x i)) x i 



< i(h*) + 3Wj] 



(15) 



Combining (13) and (15) in (12), we get 



„t+i 



w\\ >£(&*) -5W>7 



□ 



The bound on for some i now follows from Lemma [jsj Let r\ — 2(1 + \J log(l /S)/l) /y/m. Notice 

that (yi — u(w ■ Xi))xi for all i are i.i.d. 0-mean random variables with norm bounded by 1, so using Lemma 
[2J ||(1/to) — u ( w ' < Now using Lemma[5j at each iteration of algorithm GLM-tron, either 



LU 



t+1 



w\\ < \\w — w\\ 



Wr\, or s(h ) < 6Wrj. If the latter is the case, we are done. If not, since 
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||w t+1 — u;|| 2 > 0, and \\w° — u>|| 2 = ||w|| 2 < W 2 , there can be at most W 2 /(Wrf) = W/(rj) iterations before 
£{hf) < 6Wn. Overall, there is some h l such that 

In addition, we can reduce this to a high-probability bound on using Lemma |4j which is applicable 

since ||u>*|| < W. Using a union bound, we get a bound which holds simultaneously for e(/i') and e(/i*). 

B Proof of Proposition [l] 

To prove the proposition, we actually prove a more general result. Define the function class 

U = {u: [-W, W] [0, 1] : u 1-Lipschitz}. 

and 

W = {x^ (x,w) :w e R d , \\w\\ < W], 

where d is possibly infinite (for instance, if we are using kernels). 

It is easy to see that the proposition follows from the following uniform convergence guarantee: 

Theorem 3. With probability at least 1 — 5, for any fixed w £ W, if we let 

^ ra 

u = argmin — y^(u(w ■ xi) - y,) 2 , 

i=l 

and define 

^ ra 

u = argmin — > (u(w ■ x{) — Eky|xi]) 2 , 
i=i 

then 

lA |A , , _, „ / . f fdW 3 / 2 log(Wm/S)\ 1/3 IW 2 log(m/<J) 

— > \u{wxi)-u(wxi)\ <0 min ^ + \ &v ' ' 

ra f— ^ \ I \ m J V m 

To prove the theorem, we use the concept of (co-norm) covering numbers. Given a function class T 
on some domain and some e > 0, we define jV"oo(e, J-) to be the smallest size of a covering set J-' C J 7 , 
such that for any / e J. there exists some /' £ J for which sup x \f[x) — f'(x)\ < e. In addition, we 
use a more refined notion of an oo-norm covering number, which deals with an empirical sample of size m. 
Formally, define A/"oo(e, J 7 , rn) to be the smallest integer n, such that for any x\, . . . , x m , one can construct 
a covering set T' C T of size at most n, such that for any / € J, there exists some /' £ f such that 
maxj =li .., i7n \f(xi) - f'{xi)\ < e. 

Lemma 6. Assuming m, l/e,W > 1, we have the following covering number bounds: 
1- Afoo(e,U) < \2 2W /\ 

2. N 00 (e,W)<(l+ 2 -^) d . 

3. ^(e,WoW) < 2 e 2^(l + ^f) d . 

4. jVoofe.MoW.m) < f(2m+ l) 1+8W/2 / £2 



W log(m/(5) 
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Proof. We start with the first bound. Discretize [— W, W) x [0, 1] to a two-dimensional grid 
{— W + ea, efr} a =o....,2W7e.6=o,...,i/e- It is easily verified that for any function u G U, we can define a piecewise 
linear function v! ', which passes through points in the grid, and in between the points, is either constant or 
linear with slope 1, and sup,,, \u(x) — u' (x)\ < e. Moreover, all such functions are parameterized by their value 
at —W, and whether they are sloping up or constant at any grid interval afterwards. Thus, their number 
can be coarsely upper bounded as 2 2W / € /e. 

The second bound in the lemma is a well known fact - see for instance pg. 63 in [Pis99 ). 

The third bound in the lemma follows from combining the first two bounds, and using the Lipschitz 
property of u (we simply combine the two covers at an e/2 scale, which leads to a cover at scale e for U o W). 

To get the fourth bound, we note that by corollary 3 in |Zha02bj . A/"oo(e, W, m) < (2m + i)i+VK 2 /« 2 . 
Note that unlike the second bound in the lemma, this bound is dimension- free, but has worse dependence 
on W and e. Also, we have A/"oo(e,W,m) < 7Voo(e,W) < l2 2W / e by definition of covering numbers and the 
first bound in the lemma. Combining these two bounds, and using the Lipschitz property of u, we get 

-(2m+l) 1+4Vv2 / £2 2 4M/ / e . 



Upper bounding 2 AW / £ by (2r 



Lemma 7. With probability at least 1 — 8 over a sample (xi,yi) 
simultaneously for any w G W, u, v! G U, 



l) 4W/2 / e2 , the the fourth bound in the lemma follows. □ 

,(x m ,y m ) the following bounds hold 



m 

y^( u ( w ' x i) ~ Vi) 2 ~ E [( u ( w ■ x ) - y) 2 



< o 



! W 2 \og(m/8) 



^ 771 

- J2(u(w ■ Xi ) - E[y\ Xl }) 2 - E [(u(w ■ x) - E[y\x}) 2 



— \u(w ■ Xi) — u (w ■ Xi)\ — E [\u(w ■ x) — u (w ■ x) 
m 1=1 



< O 



< O 



'W 2 \og{m/8) 



'W 2 log(m/<5) 



Proof. Lemma ^ tells us that Mooi^-M o VV,m) < 2 (2m + l) 1 + 8H/ / c . It is easy to verify that the same 
covering number bound holds for the function classes {(x,y) i— >• (u(w • x) — y) 2 : u G U,w G W} and 
{x i — y (u(w ■ x) — ¥\y\x]) 2 : u G U,w G W}, by definition of the covering number and since the loss 
function is 1-Lipschitz. In a similar manner, one can show that the covering number of the function class 
{x h-> \u(w ■ x) - u'(w ■ x)\ : u,v! £U,w G W} is at most f (2m + l) 1+32W/2 / £ . 

Now, one just need to use results from the literature which provides uniform convergence bounds given a 
covering number on the function class. In particular, combining a uniform convergence bound in terms of the 
Rademacher complexity of the function class (e.g. Theorem 8 in |BM02] ). and a bound on the Rademacher 
complexity in terms of the covering number, using an entropy integral (e.g., Lemma A. 3 in [SSTlOj ). gives 
the desired result. □ 

Lemma 8. With probability at least 1 — 8 over a sample (x\, y\), . . . , (x m , y m ), the following holds simulta- 
neously for any w G W: if we let 



Uw({w, •>) = argmin — y^Mw ■ x t ) - Vif 



denote the empirical risk minimizer with that fixed w, then 



E(u w (w ■ x) - y) 2 - inf E(u(w ■ x) - y) 2 < O [W 



dlog(Wm/8)\ 2/3 
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Proof. For generic losses and function classes, standard bounds on the the excess error typically scale as 
0(1/ \Jm). However, we can utilize the fact that we are dealing with the squared loss to get better rates. In 
particular, using Theorem 4.2 in |Men02j . as well as the bound on jV tx ,(e,W) from Lemma |6j we get that for 
any fixed w, with probability at least 1 — 6, 



E(u w {w ■ x) - yf - inf E(u(w ■ x) - yf < O [W 



/log(l/<5)V /3N 



ueu \ \ m J J 

To get a statement which holds simultaneously for any w, we apply a union bound over a covering set of W. 
In particular, by Lemma|6j we know that we can cover W by a set W of size at most (1 + 2W/e) d , such that 
any element in W is at most e-far (in an oo-norm sense) from some w' € W. So applying a union bound 
over W', we get that with probability at least 1 — 6, it holds simultaneously for any w' € W that 

E(u w ,({ W ',x))-yf-ME { u(( W ',x))-yf < O (w ( ^(l/6) + dlo g (l + 2W/e) Y 3 ) . 



Now, for any w £ W, if we let w' denote the closest element in W' , then u(w ■ x) and u((w' , #)) are e-close 
uniformly for any it eM and any a;. From this, it is easy to see that we can extend (16) to hold for any W, 
with an additional O(e) element in the right hand side. In other words, with probability at least 1 — 6, it 
holds simultaneously for any w <E W that 

_ , , ... , 'log(2/ ( 5)+dlog(l + 2PF/e)^ 2 :x 

E(n u ,(w ■ x) — y) — ml E(u(iu ■ x) — y) < U \W ' 



Picking (say) e = 1/m provides the required result. □ 

Lemma 9. Let F be a convex class of functions, and let f* = arg min^gp E[(/(x) — y) 2 ]. Suppose that 
E[y|x] € U o W. Then for any f e F, it holds that 

Him - Vf] - E[(fW - Vf\ > E [{f{x) t(x)f] > (E 0/^) - f*(x)\]f . 

Proof. It is easily verified that 

E[(f(x) - yf] E[(f*(x) - y) 2 } = E x [(f(x) - E[y\x}) 2 - (f (x) - E^la;]) 2 ]. (17) 

This implies that /* = argminj 6f E[(f(x) - E[y\x}) 2 ]. 

Consider the L2 Hilbert space of square-integrable functions, with respect to the measure induced by the 
distribution on x (i.e., the inner product is defined as (/, /') = E x [f(x)f'(x)]). Note that E[y|x] € U o W is 
a member of that space. Viewing E[y|x] as a function y(x), what we need to show is that 

11/ - yf -\\r-vf> 11/ -rf. 

By expanding, it can be verified that this is equivalent to showing 

(r-yj-n >o. 



To prove this, we start by noticing that according to (IT), /* minimizes ||/ — y\\ over Therefore, for any 
/ € F and any e € (0, 1), 

||(l-e)/* + e /-2/|| 2 -||/*-y|| 2 >0, (18) 



as (1 — e)f* + ef € F by convexity of F. However, the right hand side of (18) equals 

e 2 \\f-f*f + 2e(r-y,f-r), 



so to ensure ( 18 ) is positive for any e, we must have (/* — y, f — /*) > 0. This gives us the required result, 
and establishes the first inequality in the lemma statement. The second inequality is just by convexity of 
the squared function. □ 
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Proof of Thm. [3} We bound ^ JZfcLi \u{w ■ x) — u(w ■ x)\ in two different ways, one which is dimension- 
dependent and one which is dimension independent. 

We begin with the dimension-dependent bound. For any fixed w, let u* be argmin ue ^ E{u(w ■ x) — y) 2 . 
We have from Lemma [8] that with probability at least 1 — 5, simultaneously for all w £ W, 

E(u(w ■ x) - yf E{u*{w • x) - yf < O (w f ***W m W ' 



and by Lemma [9j this implies 

E[\u(w x) - u* (w x)\] < O 



dW 3 / 2 \og(Wm/5) 



1/3N 



(19) 



Now, we note that since u* = argmin ug n E(u(w ■ x) — y) 2 , then u* = arg min ue ^ E(u(w ■ x) — E[y\x\) 2 as 
well. Again applying Lemma [8] and Lemma [9] in a similar manner, but now with respect to u, we get that 
with probability at least 1 — 5, simultaneously for all w € W, 



E[\u(wx) -u*(wx)\] < O 



Combining (19) and (20), with a union bound, we have 

E[\u{w x) - u(w ■ x)\] <0 



dW 3 / 2 \og(Wm/S) 



dW 3 / 2 \og(Wm/5) 



1/3N 



1/3N 



(20) 



Finally, we invoke the last inequality in Lemma [7J using a union bound, to get 



1 m ( 
— } \u(w ■ x) — u(w ■ x)\ < O I 

m t—f \ 

i=i \ 



fdW 3 ' 2 log(Wm/5) 
V m 



1/3 



W 2 \og{m/5) 



We now turn to the dimension-independent bound. In this case, the covering number bounds are different, 
and we do not know how to prove an analogue to Lemma [8] (with rate faster than 0(l/^/m)). This leads to 
a somewhat worse bound in terms of the dependence on m. 

As before, for any fixed w, we let u* be arg min u6 ^ E[(w(w • x) — y) 2 ) . Lemma[7]tells us that the empirical 
risk — Y^i=i( u ( w ' x i) ~ Vi) 2 is concentrated around its expectation uniformly for any u,w. In particular, 



as well as 



1 m ( 
- Y^iuiw ■ xi) - yi) 2 - E [u(w ■ x) - y) 2 } < O 

1=1 \ 

m / 
- ^2(u*(w ■ Xi) - yi) 2 - E [(u*(w ■ x) - y) 2 } < O 

i=l \ 



W 2 log(m/5) \ 
m J 

W 2 \og{m/5) 



but since u was chosen to be the empirical risk minimizer, it follows that 



E [(u(w ■ Xi ) - Vl ) 2 ] - E [(«* (wx)- y) 2 ] < O U W2 ^ m l 5 ) 



so by Lemma |9j 



E[\u* {w ■ x) - u{w x)\] <0 



W 2 log{m/6)\ 1/4 



(21) 



Now, it is not hard to see that if u* = argmin ue ^ E[(u(w-x)—y) 2 }, then u* = argmin U £w E[(u(w-x)— E[y\x]) 2 ] 
as well. Again invoking Lemma [7J and making similar arguments, it follows that 



E[\u*{w ■ x) - u(w ■ x)\] < O 



/W 2 log(m/£) 
\ m 



1/4N 



(22) 



1G 



Combining (21) and (22), we get 

E[\u(wx) -u(wx)\] < O 



W 2 log(m/<S)\ a 



We now invoke Lemma [7] to get 

1 m ( 

— y \u(w ■ — u(w ■ Xi)\ < O 

i=i \ 



fW 2 \og{m/S) 
\ m 
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